#!/usr/bin/env python
# -*- coding: iso-8859-1 -*-
from sage.all import *
#Actually, sage redefines the builtin random module with a function we don't want.
import random #@Reimport
from bits import *

def xgcd(p1, p2):
    '''
    Invokes normalized gcd on p1, and p2.
    @param p1:
    @param p2:
    '''
    d, c1, c2 = p1.xgcd(p2)
    if d.is_monic():
        return d, c1, c2
    divisor = d.monic()
    factor = d // divisor
    return divisor, c1 // factor, c2 //factor

def embeddingDegree(q, n):
    '''
    Returns the embedding degree of q, to n. That is: is the smallest positive integer k
    such that n | q^k - 1. In other words, F q^k is the smallest field
    that contains the group of the nth roots of unity.
    @param q: integer. Must be coprime to n
    @param n: integer.
    '''
    k = 1
    qExpk = q
    while (qExpk - 1) % n != 0:
        k += 1
        qExpk *= q
    return k

class Divisor(object):
    @classmethod
    def safeNew(self, u, v, curve):
        ret = Divisor(u, v, curve)
        ret.assert_isvalid()
        return ret

    def __init__(self, u, v, curve):
        self.u = u
        self.v = v
        self.curve = curve
        self.f = self.curve.f
        self.h = self.curve.h

    def toPrincipal(self):
        ret = self
        nullDivisor = self.curve.nullDivisor()
        retrec = ret.reduce()
        while not retrec == nullDivisor:
            ret += self
            retrec = ret.reduce()
        return ret
    
    def order(self):
        ret = self.reduce()
        nullDivisor = self.curve.nullDivisor()
        order = 1
        while not ret == nullDivisor:
            ret = ret.addReduce(self)
            order += 1
        return order
    
    def addReduce(self, o):
        return (self + o).reduce()

    def div(self, u, v):
        return Divisor(u, v, self.curve)


    def __hash__(self):
        return hash(self.asTuple())


    def __eq__(self, o):
        return self.u == o.u and self.v == o.v
    
    def __ne__(self, o):
        return not (self == o) 

    def __rmul__(self, integer):
        return self.times(integer)


    def times(self, i):
        '''
        Returns i * self (integer pre-multiplication of a group element)
        @param i: an integer
        '''
        if i == 0:
            return self.curve.nullDivisor()
        if i < 0:
            return self._times_(-i)
        return self._times_(i)

    def _times_(self, i):
        ret = self
        for Si in bitsWithoutLeadingBit(i):
            ret += ret
            if Si == 1:
                ret += self
        return ret


    def assert_isvalid(self):
        '''
        Assertions that this divisor represents u divisor on the picard group 
        of the hyperellptic curve defined by self.curve
        '''
        genus = self.curve.genus()
        u, v = self
        f, h = self.f, self.h
        assert u.degree() <= genus, "bad degrees: %d is greater than and %d" % (u.degree(), genus)
        assert u.divides(v ** 2 + v * h - f), "%r does not divide %r" % (u, v ** 2 + v * h - f)
        assert v.degree() < u.degree(), "bad degrees: u and v"
        return True
    
    

    def asTuple(self):
        return self.u, self.v

    def __iter__(self):
        return iter(self.asTuple())

    def __add__(self, divisor):
        u1, v1 = self
        u2, v2 = divisor
        d1, e1, e2 = xgcd(u1, u2)
        d, c1, c2 = xgcd(d1, v1 + v2 + self.h)
        s1 = c1 * e1
        s2 = c1 * e2
        s3 = c2
        u = (u1 * u2) // (d ** 2)
        v = ((s1 * u1 * v2 + s2 * u2 * v1 + s3 * (v1 * v2 + self.f)) // d).mod(u)
        return self.div(u, v)

    def __getitem__(self, index):
        return self.asTuple()[index]

    def __repr__(self):
        return "[%r, %r]" % self.asTuple()

    def __neg__(self):
        return self.div(self.u, -self.h - self.v)

    def __sub__(self, divisor):
        return self + (-divisor)

    def reduce(self):
        u, v = self
        tu = (self.f - v * self.h - v ** 2) // u
        tv = (-self.h - v).mod(tu)
        while tu.degree() > self.curve.genus():
            u = tu
            v = tv
            tu = (self.f - v * self.h - v ** 2) // u
            tv = (-self.h - v).mod(tu)
        return self.div(tu.monic(), tv)

class Curve(object):
    '''
    Implements algorithms for on Hyperlliptic Curve Cryptography
    '''

    @classmethod
    def createSimpleCurve(self, f, polynomialRing):
        return Curve(f, polynomialRing(0), polynomialRing)

    def __init__(self, f, h, polynomialRing):
        self.f = f
        self.h = h
        self.polynomialRing = polynomialRing
        self.curve = HyperellipticCurve(f, h)
        self.x = polynomialRing.gen()
        self.field = polynomialRing.base_ring()
        self.one = polynomialRing(1)
        self.zero = polynomialRing(0)
        self.fieldOne = self.field(1)

    def _div(self, x, y):
        '''Returns a div in Mumford representation. Does not check if [x,y] is valid'''
        return Divisor(x, y, self)


    def div(self, x, y):
        '''Returns a div in Mumford representation'''
        return Divisor.safeNew(x, y, self)

    def divisorOf(self, arg):
        '''Returns a div in Mumford representation from a list or tuple'''
        u, v = arg
        return self.div(u, v)

    def nullDivisor(self):
        return Divisor(self.one, self.zero, self)

    def genus(self):
        return (self.f.degree() - 1) / 2

    def asPolynomial(self, element):
        return self.polynomialRing(element)

    def fhValuesOf(self, element):
        return (self.f(element), self.h(element))

    def randomDivisor(self):
        '''
        Implementation of Random Divisors based on the one given by Niels Lubbes thesis
        (page 78 of by http://www.risc.uni-linz.ac.at/people/nlubbes/mt/mt/hcdlpV2.pdf)
        '''
        degree = self.genus()
        D = self.nullDivisor()
        hashtable = {}
        for i in xrange(degree):
            y = None
            element = None
            while y is None:
                element = self.field.random_element()
                if element not in hashtable:
                    beta, alfa = self.fhValuesOf(element)
                    pol = self.x ** 2 + alfa * self.x - beta
                    roots = pol.roots(multiplicities = False)
                    if len(roots) > 0:
                        y = roots[0]
                        hashtable[element] = y
                else:
                    y = hashtable[element]
            Dp = self._div(self.x - element, self.asPolynomial(y))
            D += Dp
            D.reduce()
        r = random.randint(1, 128)
        return (r * D).reduce()

    def randomPrincipalDivisor(self):
        return self.randomDivisor().toPrincipal()

    
    def miller(self, S, d, D1, D2):
        try:
            return self.internaMiller(S, d, D1, D2)
        except (ZeroDivisionError, NonDisjointSupport):
            return self.fallbackMiller(S, d, D1, D2)
    
    def fallbackMiller(self, S, d, D1, D2):
        while True:
            r = self.randomDivisor()
#            print "the random divisor is %r" % r
            Da = (D2 + r).reduce()
            Db = r
            try:
                second = self.internaMiller(S, d, D1, Da)
                first = self.internaMiller(S, d, D1, Db)
                return first / second
            except (ZeroDivisionError, NonDisjointSupport):
                pass
    
    def shifted_miller(self, S, d, D1, D2, r):
        '''
        Debugging vesion of internaMiller non effective, that allows the use of a specific shift
        '''
        print "the random divisor is %r" % r
        Da = (D1 + r).reduce()
        Db = r
        return self.internaMiller(S, d, Da, D2) / self.internaMiller(S, d, Db, D2)



    def internaMiller(self, S, d, D1, D2):
        '''
        Computes f_{s,D1}(D2). Inspired on the paper Hyperellpitic Pairings:
        http://www.isg.rhul.ac.uk/~sdg/Pairing2007-survey.pdf
        
        @param S: Order of D1.
        @param d: Final exponent. On a finite field of characteristic q,
            this must be (q ** embeddingDegree(q, S) - 1)/S
        @param D1: divisor of the curve. See methods div and divisorOf.
        S * D1 must be nullDivisor.
        @param D2: divisor of the curve 
        '''
        D = D1
        u2, v2 = D2 #@UnusedVariable: v2
        f1, f2, f3 = self.one, self.one, self.fieldOne
        for Si in bitsWithoutLeadingBit(S):
            f1 = (f1 * f1).mod(u2)
            f2 = (f2 * f2).mod(u2)
            f3 = f3 * f3
            D, (h1, h2, h3) = self.millerstep(D, D, D2)
            f1 = (f1 * h1).mod(u2)
            f2 = (f2 * h2).mod(u2)
            f3 = f3 * h3
            if Si == 1:
                D, (h1, h2, h3) = self.millerstep(D, D1, D2)
                f1 = (f1 * h1).mod(u2)
                f2 = (f2 * h2).mod(u2)
                f3 = f3 * h3
        f0 = u2.resultant(f1) / (f3 ** (u2.degree()) * u2.resultant(f2))
        return self.assertNonZero(f0 ** d)
    
    def assertNonZero(self, value):
        if value == self.zero:
            raise NonDisjointSupport()
        return value

    
    def plotFile(self, filename, xRange = 2, yRange = 3, precisao = 350):
        x = polygen(QQ) #@UnusedVariable, used on eval
        f = eval(str(self.f).replace('^', '**')) + 0*x
        h = eval(str(self.h).replace('^', '**')) + 0*x
        func = lambda px, py: py ** 2 + h(px) * py - f(px)
        c = contour_plot(func, (-xRange, xRange), (-yRange, yRange),
                         plot_points = precisao, cmap = 'winter', fill = False, contours = [0])
        R = line([(0, -yRange), (0, yRange)], linestyle = '--', rgbcolor = (0, 0, 0))
        L = line([(-xRange, 0), (xRange, 0)], linestyle = '--', rgbcolor = (0, 0, 0))
        graph =  c + L + R
        graph.save(filename)
        


    def millerstep(self, d1, d2, e):
        u1, v1 = d1
        u2, v2 = d2
        ue, ve = e
        d1, e1, e2 = xgcd(u1, u2)
        d, c1, c2 = xgcd(d1, v1 + v2 + self.h)
        h1til = d.mod(ue)
        h2til = self.one
        h3 = self.fieldOne
        s1 = c1 * e1
        s2 = c1 * e2
        s3 = c2
        u = (u1 * u2) // (d ** 2)
        v = ((s1 * u1 * v2 + s2 * u2 * v1 + s3 * (v1 * v2 + self.f)) // d).mod(u)
        while u.degree() > self.genus():
            ulinha = ((self.f - v * self.h - v * v) // u).monic()
            vlinha = (-self.h - v).mod(ulinha)
            h1til = (h1til * (ve - v)).mod(ue)
            h2til = (h2til * ulinha).mod(ue)
            if v.degree() > self.genus():
                h3 = -(v.leading_coefficient()) * h3
            u = ulinha
            v = vlinha
        return ((u, v), (h1til, h2til, h3))

    def __str__(self):
        return "Curve = : %s" % self.curve


class NonDisjointSupport(Exception):
    pass